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Abstract 

The computational fluid dynamics code FLUENT was used to study Rayleigh 
instability at large temperature differences in a sealed gas-filled enclosure with a cold top 
surface and a heated bottom wall (Benard problem). Both steady state and transient 
calculations were performed. The results define the boundaries of instability in a system 
depending on the geometry, temperature and pressure. It is shown that regardless of how 
fast the bottom-wall temperature can be ramped up to minimize the time spent in the 
unstable region of fluid motion, the eventual stability of the system depends on the 
prevailing final pressure after steady state has been reached. Calculations also show that 
the final state of the system can be different depending on whether the result is obtained via 
a steady-state solution or is reached by transient calculations. Changes in the slope of the 
pressure-versus-time curve are found to be a very good indicator of changes in the flow 
patterns in the system. 
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Nomenclature 

C p heat capacity at constant pressure (kJ/kg K) 

g gravitational constant (m/s 2 ) 

Gr Grashof number (dimensionless) 

k thermal conductivity (W/m K) 

L characteristic length (m) 

P pressure (atm) 

Pi initial pressure (atm) 

Pf final pressure (atm) 

Pr Prandtl number (dimensionless) 

Q total heat transfer rate (J/s) 

Ra Rayleigh number (dimensionless) 

Ra c critical Rayleigh number (dimensionless) 

T temperature (K) 

T average temperature (K) 

T c cold-wall temperature (K) 

Th hot- wall temperature (K) 

T; initial temperature (K) 

V max maximum velocity (m/s) 

Greeks 

a thermal diffusivity (m 2 /s) 

[3 t coefficient of thermal expansion at constant pressure (K 1 ) 
p dynamic viscosity (kg/m s) 

v kinematic viscosity (m 2 /s) 

p density (kg/m 3 ) 

p t initial density (kg/m 3 ) 
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1. Introduction 


When a horizontal enclosed space is heated to a uniform bottom- wall temperature 
greater than that of the top wall, an upwardly decreasing temperature profile is established. 
For a sufficiently small temperature difference between the bottom and top walls, viscous 
forces are greater than buoyancy forces, and the fluid remains quiescent. However, as this 
temperature difference increases, a point will be reached at which the buoyancy forces will 
become greater than the viscous forces, and the fluid begins to move. This is known as the 
classic Benard problem. 

The Rayleigh number is used to define the onset and development of buoyant flow 
under these conditions. It is defined as: 

Ras (p,- P ) ge (1) 

Pi av 

where 



( 2 ) 


The Rayleigh number can also be expressed as the product of the Grashof number, Gr, and 
the Prandtl number, Pr: 

Ra = Grx Pr (3) 

where 


Gr = 


(a ~p)g^ 

Pi v 1 


(4) 


and 
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( 5 ) 


Pr s = p£? 

a k 

The fluid remains stable if the Rayleigh number is less than the critical Rayleigh number, 
Ra c . Once the critical Rayleigh number is exceeded, finger-like intrusions of the lighter 
gas into the heavier gas and the heavier gas into the lighter gas result. Goldstein and 
Volvino (1995) have presented a summary of various investigations in the literature for the 
onset of convection and flow development in horizontal layers. Widely ranging values of 
the critical Rayleigh number were reported, anywhere from 1000 to 3000. 

While the above studies used different criteria and different experimental 
observations as the basis for determining the onset of instability, one possible explanation 
for this wide range of values is that the Rayleigh and Grashof numbers are derived 
assuming the Boussinesq approximation. This assumes that all fluid properties are constant, 
except for density in the buoyancy term of the momentum equation, and that viscous 
dissipation can be neglected (Kays and Crawford, 1980). Under these conditions, the 
density difference can be approximated as: 


(Pt ~P) 


= Pr(T-Ti) 


where 


P) V8T 


which is the coefficient of thermal expansion at constant pressure. 

While the temperature dependence of thermodynamic and transport properties can 
often be safely neglected at small temperature differences, they become important at large 
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temperature differences and their effects must be included in order to obtain accurate 
descriptions of the state of flow inside the enclosure. This is particularly true for gases, 
which tend to become more stable at higher temperature differences due to the dampening 
effect of gas viscosity at high temperatures. Under non-Boussinesq conditions, the 
predictive abilities of dimensionless numbers derived by assuming Boussinesq conditions 
require more attention. Gray and Giorgini (1976) and Paolucci and Chenoweth (1987) 
reported on the validity of the Boussinesq approximation and departures from it. These 
studies led to the specification of conditions under which it is justified to neglect certain 
terms in the Boussinesq approximation. 

In a previous paper, Gokoglu and Kuczmarski (2003) showed via a steady state 
analysis that for temperature differences larger than justifiable for the Boussinesq 
approximation, there are situations were a gaseous system will never experience the onset 
of buoyant convection, no matter how large the temperature difference becomes at a given 
pressure. This was due to the temperature-dependent gas properties, particularly that of the 
kinematic viscosity, which tended to counteract the tendency toward instability at higher 
temperatures. In addition, two critical hot surface temperatures were shown to exist for a 
given pressure, above and below which the system was stable. 

In this paper, transient calculations have been performed to examine the behavior of 
the system as it pressurizes while the bottom wall is heated at different rates. These 
transient solutions are carried out over a sufficiently long time period to allow comparison 
between results obtained from transient solutions and those obtained using steady-state 
solutions at the final prevailing pressure. The development of the flow field over time for 
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different heating rates and its effect on the temperature field is demonstrated. In addition, 
the effect of varying the enclosure height on the predicted range of instability is examined. 

2. Numerical Models 

This investigation used the computational fluid dynamics code FLUENT, which 
uses a finite volume method to discretize the continuity, momentum, and energy equations. 
Two 2-D models of a rectangular geometry were constructed, as shown in Fig. 1. Case 2 
has a distance between the hot and cold walls that is twice that of case 1 . Adiabatic 
boundary conditions were set along both side walls, and the top cold-surface temperature 
was set at 300 K. The bottom hot-surface temperature was set at 2500 K for the steady 
state cases, or ramped to a final temperature of 2500 K for the transient cases. Neon gas 
was used in the enclosure. For the steady-state cases, the operating pressure was defined, 
along with a hot-wall temperature of 2500 K. The standard discretization scheme in 
FLUENT was used for pressure and the SIMPLE algorithm for the pressure-velocity 
coupling. A second order upwind scheme was used for the continuity, momentum, and 
energy equations. For the transient calculations, an initial pressure was specified, but a 
floating operating pressure was defined in the model to allow the pressure to increase as the 
gas was heated. A second order upwind scheme was used for the momentum and energy 
equations, and for the density interpolation scheme. The PRESTO! discretization scheme 
was used for pressure, and the PISO algorithm was employed for the pressure-velocity 
coupling. A second order upwind scheme was used for the continuity, momentum, and 
energy equations. A time step of lxlO" 4 second was used for all transient calculations. 
Larger time steps led to problems with convergence. The under-relaxation factors were set 
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to 1 for pressure, density, body forces, momentum, and energy. This helped speed the 
computations that were slowed by the necessity of choosing a small time step. The 
solutions were considered converged for steady state solutions when the scaled residuals 
dropped below lxlO" 3 for the continuity equation, lxlO" 5 for the momentum equations, and 
lxlO" 6 for the energy equation. Convergence for the transient solutions required these 
values to be obtained for each time step. All transient cases were run up to ten seconds. 

Four different grid densities were examined for case 1 to study the solution grid 
sensitivity. The total heat transfer rate at the heated bottom wall was compared among the 
four grids, since this heat transfer into the system was the ultimate source for all changes 
experienced by the system. Table 1 shows the number of computational cells used for each 
grid, along with the total heat transfer rate at the bottom wall, Q, and the change in this rate 
between consecutive grids. The change in the total heat transfer rate is of the same order of 
magnitude between grids 1 and 2 as it is between grids 2 and 3; both are at very small 
values. Because the computational time increases with the grid density, grid 2 was chosen 
to balance higher resolution with reasonable run times. All results in the rest of this paper 
for case 1 were obtained using grid 2. 

For case 2, the grid used was based on grid 2 discussed above. The number of 
computational cells in the X-direction remained the same at 240 cells, but the number of 
computational cells in the Y-direction was simply doubled from 40 cells to 80 cells since 
the height of case 2 is double that of case 1 . 
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3. Results 


First, the upper and lower temperature limits on stability were determined for the 
two cases with different enclosure heights by running steady-state cases at various pressure 
levels. The results are shown in Fig. 2. It can be seen that regardless of the temperature 
difference, some systems always remain stable below a certain minimum pressure, 
depending on the height of the geometry. Neon was stable for all hot-surface 
temperatures at 1 atm of pressure or less for case 1, but the pressure needed to be decreased 
to 0.36 atm or less in case 2 before neon was stable at all hot-surface temperatures. In 
addition, the range of pressures and temperatures for which the system is unstable is 
expectedly much greater for case 2 than for case 1 since the characteristic length, i.e., the 
height between the top and bottom walls in this study, is doubled in case 2, resulting in an 
eight-fold increase in the Rayleigh number. These results also confirm the original findings 
of Gokoglu and Kuczmarski (2003) that the tendency of a system to become more unstable 
can actually decrease when hot-surface temperature is increased beyond a certain level 
where the dampening effect of gas viscosity becomes more prevalent. 

Next, transient runs were conducted to examine the behavior of the system as it 
dynamically develops. The same scales were used for Figs. 8-17, i.e., a range of 2xl0~ 8 to 
5x10"' m/s for velocity and a range of 300 to 2500K for temperature. The curves showing 
the rise in pressure with the bottom-wall temperature are plotted in Figs. 3 and 4 for both 
cases at initial pressures of 1.0 and 1.5 atm. The bottom-wall temperature was ramped to a 
final temperature of 2500 K over 0.1, 0.5, 1.0, and 5.0 seconds. The curves are 
superimposed over the results obtained from steady-state calculations showing the stable 
and unstable regions. For all ramp times, the system crosses over into the unstable region 
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at least part of the time. For case 1, at an initial pressure of 1 atm, all cases cross back over 
to the stable region at the end of the ramp, although they remain very close to the transition 
line. For case 2, at an initial pressure of 1 atm, the system quickly crosses into the unstable 
region and remains there for the 0.5, 1.0, and 5.0 second ramp times. For the 0.1 second 
ramp time, the curve crosses over into the stable region again, but after the final hot-wall 
temperature is reached, the pressure continues to rise towards its steady state, and the 
system ends up back in the unstable region. So, even though the bottom-wall temperature 
can be ramped at a rate that causes the system to undergo a stable-unstable- stable transition, 
the system will eventually end up in a region defined by the final pressure the system rises 
to. This pressure can be estimated using the ideal gas law: 



nRT 

V 


(8) 


The number of moles of gas in the system can be determined from the initial conditions, but 
the temperature to use in this equation is not as obvious. The arithmetic mean of the hot- 
and cold-wall temperatures in the system, customarily used to evaluate gas properties, can 
be used (Kays and Crawford, 1980; Bird, Stewart, and Lightfoot, 1963): 


T 


T +T 
2 


(9) 


Table 2 shows the final system pressure calculated from equation (8), the pressure predicted 
by FLUENT at 10 seconds, and the percent differences between the two, determined by the 
equation: 
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The largest difference observed is less than 15%, a relatively close agreement. 

The pressure rise with time for various bottom-wall temperature ramp rates is 
shown in Figs. 5 and 6. A change in the slope of the curve is a good indicator of a change 
in the flow field. This is most clearly illustrated in Fig. 6b for case 2, with an initial 
pressure of 1.5 atm and a bottom- wall temperature ramp rate of 5.0 sec. At around the five 
second mark, the pressure curve changes from a steady rise to a brief plateau. Fig. 7a 
shows this region in more detail with higher resolution. During about a 0.6-second interval 
where the pressure drops, Fig. 7b - 7e shows the transition in the velocity field from six 
convective cells down to four. An accompanying change in heat transfer rates and 
temperature in the enclosure causes the observed change in the pressure curve. This 
observation can be a convenient experimental tool for diagnosing the expected changes in 
flow behavior in such systems. 

Regarding case 1, it was shown in Fig. 3(a) that the pressure vs. hot- wall 
temperature curves end up in the stable region near the transition between unstable and 
stable flow for all hot-wall temperature ramp times investigated. This stability map was 
obtained using steady-state calculations. Fig. 8 shows that when a transient solution is 
obtained, convective rolls are predicted to be present in the system for the 0.1 and 0.5 sec. 
ramp times. However, the velocity is at such low values at selected times during the 
transient run that the temperature fields are not significantly affected. This is in contrast to 
the results obtained from a transient solution shown in Fig. 9. Convective rolls are 
predicted to be present at a sufficient strength to affect the temperature fields. In both 
figures, the transient solution is seen to no longer be changing at ten seconds, indicating 
that a steady state has been reached. It can be seen that the transient results at ten seconds 
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are different than those predicted by the steady state solution, pointing out that a transient 
solution may be necessary to accurately predict the final state of a system such as this. 

Figs. 10 and 11 show similar results between the steady-state analysis and the transient 
calculations carried out to ten seconds. Compared to the cases with an initial pressure of 1 
atm, these cases with the initial pressure of 1.5 atm show a much greater effect on the 
temperature field from the velocity field. 

Regarding case 2, Figs. 12, 13, and 14 show that the greater height in case 2 has led 
to higher velocities compared to case 1, with a corresponding effect in the temperature 
field. Similar results for an initial pressure of 1.5 atm are shown in Figs. 15, 16, and 17. At 
an initial pressure of 1 atm, the number and/or size of the predicted convective rolls is 
different between the steady-state solution and the transient solution taken out to ten 
seconds. For the cases with an initial pressure of 1.5 atm, on the other hand, steady-state 
calculations can sometimes predict the final state of the system consistently with the 
transient calculations. Such discrepancies between steady-state and transient calculations 
points out again that the history of the development of flow structures and their 
corresponding temperature fields have a profound effect on the eventual final state which 
can not be properly captured by only steady-state calculations. 

4. Conclusions 

This study demonstrated a capability to define the range of instability in a gaseous 
enclosed system depending on its geometry, temperature, and pressure. It is shown by 
transient calculations that while the bottom-wall temperature can be ramped fast enough to 
minimize the time spent in the unstable region of fluid motion, the final state of the system 
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depends on the final pressure eventually attained after steady state has been reached. The 
results also show that steady-state and transient solutions can give different results as to the 
final state of the system indicating that the history of the development of flow instabilities 
is important and that steady-state calculations can not always be relied upon. It is also 
demonstrated that changes in the monotonic behavior of the pressure versus time curve in a 
transient treatment of the problem is a good indicator of changes in flow patterns in the 
system and can also be used as an experimental diagnostic tool. 

Many applications involving heat transfer in gaseous media, such as chemical vapor 
deposition, combustion, and cryogenics, require knowledge of the onset of convective 
instability in order to design systems to perform as needed. This work has shown that 
high temperatures can be effectively used to actually suppress buoyant convection in 
gaseous systems utilizing the enhanced dampening effect of increased gas viscosities. 

Using existing heat transfer correlations involving the Rayleigh and Grashof numbers under 
conditions where the Boussinesq approximation is no longer valid, a common practice, 
could be seriously misleading since the system stability can not correctly be predicted by 
these dimensionless numbers. In addition, the intensity of the buoyant convection can be 
misjudged by not considering the history of the system during heating. Depending on the 
system conditions, the temperature dependent properties of the gas must be included to 
fully understand the heat transfer occurring in the system. 
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Fig. 1: Diagrams of configurations modeled for this paper, (a) case 1, (b) case 2, which 
is the same in length as case 1 , but double the height. 
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Fig.2: a) A schematic of calculated stable and unstable regions for neon for (a) case 1 
(b) case 2. Hot surface temperatures at which transition occurs at various pressures 
are indicated. 
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Fig.3: Case 1: pressure vs. hot wall temperature for various hot wall temperature ramp 
times, (a) P ± = 1.0 atm, (b) P ; = 1.5 atm. 
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Fig.4: Case 2: pressure vs. hot wall temperature for various hot wall temperature ramp 
times, (a) P ± = 1.0 atm, (b) P ; = 1.5 atm. 
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Fig.5: Case 1: pressure vs. time for various hot wall temperature ramp times, 
(a) Pj = 1.0 atm, (b) P ; = 1.5 atm. Note the different scales on the two plots. 
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Fig.6: Case 2: pressure vs. time for various hot wall temperature ramp times, 
(a) Pj = 1.0 atm, (b) Pj = 1.5 atm. Note the different scales on the two plots. 





time (s) 

(a) 






Fig.7: Case 2: (a) pressure vs. time for the five second temperature ramp between 4.5 and 
6.5 s; velocity fields at: (b) 5.0 s, ( c) 5.2 s, (d) 5.4 s, (e) 5.6 s. 
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Fig. 8: Velocity and temperature fields in neon, case 1. Final T h = 2500 K, P, = 1.0 atm. 
Transient solution with a 0.1 second bottom wall temperature ramp after: 

(a) 0.3 s (V max =0.007 m/s), (b) 0.5 s (V max = 0.004 m/s), (c) 2.0 s (V max = 0.0007 m/s), 
(d) 5.0 s (V max = 0.0006 m/s), (e) 10.0 s (V max = 0.0006 m/s).; (f) steady state solution 
(V m ax = 0-009 m/s). 0.5 second temperature ramp yields similar results. 
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Fig. 9: Velocity and temperature fields in neon, case 1. Final T h = 2500 K, P, = 1.0 atm. 
Transient solution with a 1 .0 second bottom wall temperature ramp after: 

(a) 0.5 s (V max =0.01 m/s), (b) 0.8 s (V max = 0.02 m/s), (c) 1.5 s (V max = 0.04 m/s), 

(d) 5.0 s (V max = 0.1 m/s), (e) 10.0 s (V max = 0.1 m/s).; (f) steady state solution 
(V m ax = 0.0.009 m/s). 5.0 second temperature ramp yields similar results. 
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Fig. 10: Velocity and temperature fields in neon, case 1. Final T h = 2500 K, P, = 1.5 atm. 
Transient solution with a 0.1 second bottom wall temperature ramp after: 

(a) 0.2 s (V max =0.03 m/s), (b) 0.4 s (V max = 0.09 m/s), (c) 0.6 s (V max = 0.2 m/s), 

(d) 0.8 s (V max = 0.2 m/s), (e) 10.0 s (V max = 0.2 m/s).; (f) steady state solution 
(V m ax = 0.2 m/s). 0.5 and 1.0 second temperature ramps yield similar results. 
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Fig. 1 1 : Velocity and temperature fields in neon, case 1 . Final T h = 2500 K, P, = 1 .5 atm. 
Transient solution with a 5.0 second bottom wall temperature ramp after: 

(a) 0.6 s (V max =0.004 m/s), (b) 0.8 s (V max = 0.01 m/s), (c) 1.1 s (V max = 0.1 m/s), 

(d) 1.3 s (V max = 0.2 m/s), (e) 10.0 s (V max = 0.2 m/s).; (f) steady state solution 

( V max = 0-2 m/s). 





Fig. 12: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.0 atm. 
Transient solution with a 0.1 second bottom wall temperature ramp after: 

(a) 0.2 s (V max =0.06 m/s), (b) 0.3 s (V max = 0.2 m/s), (c) 0.5 s (V max = 0.3 m/s), 

(d) 0.8 s (V max = 0.3 m/s), (e) 10.0 s (V max = 0.3 m/s).; (f) steady state solution 
(V m ax = 0-3 m/s). 0.5 second temperature ramp yields similar results. 




Fig. 13: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.0 atm. 
Transient solution with a 1.0 second bottom wall temperature ramp after: 

(a) 0.4 s (V max =0.03 m/s), (b) 0.6 s (V max = 0.2 m/s), (c) 0.8 s (V max = 0.3 m/s), 

(d) 3.0 s (V max = 0.3 m/s), (e) 10.0 s (V max = 0.3 m/s).; (f) steady state solution 
(V max = 0.3 m/s). 




Fig. 14: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.0 atm. 
Transient solution with a 5.0 second bottom wall temperature ramp after: 

(a) 0.6 s (V max =0.006 m/s), (b) 1.0 s (V max = 0.1 m/s), (c) 1.5 s (V max = 0.2 m/s), 

(d) 2.0 s (V max = 0.3 m/s), (e) 10.0 s (V max = 0.3 m/s).; (f) steady state solution 
(V max = 0.3 m/s). 




Fig. 15: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.5 atm. 
Transient solution with a 0.1 second bottom wall temperature ramp after: 

(a) 0.2 s (V max =0.07 m/s), (b) 0.4 s (V max = 0.5 m/s), (c) 0.5 s (V max = 0.4 m/s), (d) 0.6 s 
(V ma x = 0-4 m/s), (e) 10.0 s (V max = 0.4 m/s).; (f) steady state solution (V max = 0.4 m/s). 
0.5 second temperature ramp yields similar results. 



Fig. 16: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.5 atm. 
Transient solution with a 1.0 second bottom wall temperature ramp after: 

(a) 0.4 s (V max =0.04 m/s), (b) 0.6 s (V max = 0.4 m/s), (c) 0.7 s (V max = 0.3 m/s), (d) 0.8 s 
(V m ax = 0-4 m/s), (e) 10.0 s (V max = 0.4 m/s).; (f) steady state solution (V max = 04 m/s). 
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Fig. 17: Velocity and temperature fields in neon, case 2. Final T h = 2500 K, P, = 1.5 atm. 
Transient solution with a 5.0 second bottom wall temperature ramp after: 

(a) 0.6 s (V max =0.008 m/s),(b) 0.9 s (V max = 0.09 m/s), (c) 3.0 s (V max = 0.3 m/s), (d) 6.0 s 
(V m ax = 0-3 m/s), (e) 10.0 s (V max = 0.3 m/s).; (f) steady state solution (V max = 0.3 m/s). 


Table 1 

Grid Sensitivity Study Results 


grid 

# cells in X 

# cells in Y 

% change in Q (J/s) 

grid 0 

60 

10 


grid 1 

120 

20 

(Q1-Q0)/Q0 = 0.03 

grid 2 

240 

40 

(Q2-Q1J/Q1 = 
0.009 


grid 3 


480 


80 


(Q3-Q2)/Q2 = 
0.002 























Table 2 

Comparison of Final Pressure Predicted by the Ideal Gas Law and FLUENT 


Case 

Pi 

(atm) 

Ramp 

time 

(sec) 

P f (atm.) 
Ideal gas 
law 

P f (atm.) 

FLUENT 

% 

difference 

1 

1.0 

0.1 

4.67 

4.27 

- 8.6 

1 

1.0 

0.5 

4.67 

4.27 

- 8.6 

1 

1.0 

1.0 

4.67 

4.17 

-10.7 

1 

1.0 

5.0 

4.67 

4.17 

-10.7 

1 

1.5 

0.1 

7.0 

6.20 

-11.4 

1 

1.5 

0.5 

7.0 

6.19 

-11.6 

1 

1.5 

1.0 

7.0 

6.19 

-11.6 

1 

1.5 

5.0 

7.0 

6.19 

-11.6 

2 

1.0 

0.1 

4.67 

4.45 

- 4.7 

2 

1.0 

0.5 

4.67 

4.45 

- 4.7 

2 

1.0 

1.0 

4.67 

3.99 

-14.6 

2 

1.0 

5.0 

4.67 

4.13 

-11.6 

2 

1.5 

0.1 

7.0 

6.81 

- 2.7 

2 

1.5 

0.5 

7.0 

6.81 

- 2.7 

2 

1.5 

1.0 

7.0 

6.81 

- 2.7 

2 

1.5 

5.0 

7.0 

6.23 

-11.0 







































































































